In this homework you will implement one of the most basic linear solvers. Gaussian elimination with bakwards substitution with row swap, you can find the details in Alg 6.1 in your textbook (or https://en.wikipedia.org/wiki/Gaussian_elimination)
Q1.a You will write the Gaussian elimination algorithm. We will proceed by writing the function rowSwap!. The inputs of the function will be:
The function will not have any output. The swaping has to be done in the matrix M. Finally, your function will raise an error if $i$ or $j$ are out of range.
In Julia, all the input variariables are passed by reference. Then you can modify the input variables outside the scope of the function. By convention, each time that a function modifies a variable outside its scope the function contains an exclamation mark.
In [7]:
function rowSwap!(A, i,j)
# input: A matrix
# i,j the row indexes to be swapped
n = size(A)[1]
# checking that i and j are within the permitted range
(i > n || j > n ) && error("Index out of range in the rowSwap")
# if the index are not the same
if i != j
buffer = A[i,:]
A[i,:] = A[j,:]
A[j,:] = buffer
end
end
Out[7]:
Q1.b You will write a function that performs the Gaussian elimination. The function takes as input:
Your function will create the augmented system, and perform the Gaussian elimination.
The output of the function will be the tuple (U,b1).
U is the upper triangular matrix resulting from the elimination and b1, is the resulting vector.
To obtain $U$ and $b1$ from your augmented matrix you perform a slicing (i.e. use [:,1:n]).
You will use your rowSwap function acting on the augmented matrix.
Finally, your function will raise an error if:
Hints:
In [21]:
function gaussianElimination(A,b)
#input: A squared matrix
# b a vector
#output: U upper triangular matrix
# b1 the resulting vector
# safety checks
(n,m) = size(A)
(n != m) && error("Matrix is not square \n")
(n != size(b)[1]) && error("Dimension mismatch \n")
# create the augmented matrix
M = hcat(A,b)
for i = 1:(n-1)
# look for the first non zero entry in the ith column
# find the indices of the non zero elements (up to machine precision)
indexj = find(abs(M[i:end,i]).> eps(1.0))
# if indexj is empty then the matrix is singular and raise error
isempty(indexj) && error("The matrix is singular \n")
# call row swap
rowSwap!(M,i,(i-1)+indexj[1])
# for loop for eliminating unknows
for j = i+1:n
M[j,:] = M[j,:] - (M[j,i]/M[i,i])*M[i,:]
end
end
# checking the last pivot!!
abs(M[n,n]) < eps(1.0) && error("The matrix is singular \n")
# slicing the matrices
U = M[:,1:n]
b1 = M[:,n+1:end]
return (U,b1)
end
Out[21]:
Q1.c You will write a function that performs the backward substitution.
The input of your function will be:
The output will be the solution to the system $Ux = b$.
Your function needs to have safeguards against a size mismatch (i.e., the size of the matrix and your vector are not compatible, or your matrix is not squared).
Hint: if you need to run a foor loop that goes from n-1 to 1, you can use the syntax
for j = n-1:-1:1
In [9]:
function backwardSubstitution(U,b)
# input: U upper triangular matrix
# b vector
# output: x = U\b
# checks for sizes
(n,m) = size(U)
(n != m) && error("Upper triangular matrix is not square \n")
(n != size(b)[1]) && error("Dimension mismatch \n")
# creating x and running the backward substitution
x = zeros(b)
x[n] = b[n]/U[n,n]
for i = (n-1):-1:1
x[i] = (b[i] - dot(U[i,i+1:end][:],x[i+1:end]))/U[i,i]
end
# returning x
return x
end
Out[9]:
You can test that your function is correct by running the following script:
In [10]:
# size of the Matrix
m = 100
# creating an upper triangular Matrix
U = diagm(m*ones(m,1)[:], 0) + diagm(rand(m-1,1)[:], 1) + diagm(rand(m-2,1)[:], 2)
# creating the rhs
b = rand(size(U)[1],1)
@time x = backwardSubstitution(U,b)
print("Residual of the backward substitution is ", norm(U*x -b)/norm(b),"\n")
The residual should be extremely small (around epsilon machine)
In [11]:
function backwardSubstitutionLoop(U,b)
# checks for sizes of U and b
n = size(U)[1]
# first step
x = zeros(b)
x[n] = b[n]/U[n,n]
for i = (n-1):-1:1
dummy = 0
for j = i+1:n
dummy += U[i,j]*x[j]
end
x[i] = (b[i] - dummy)/U[i,i]
end
return x
end
Out[11]:
In [25]:
function solveGauss(A,b)
# input: A square matrix
# b vector
# output: x = A\b
(U,b1) = gaussianElimination(A,b)
return backwardSubstitution(U,b1)
end
Out[25]:
You can test your code by running the following script
In [26]:
# size of the Matrix
m = 100
# creating the Matrix
A = randn(m,m) + m*eye(m)
println("The conditioning number of A is " ,cond(A))
# creating the rhs
b = rand(size(A)[1],1)
@time x = solveGauss(A,b)
print("Residual of the solver is ", norm(A*x -b)/norm(b),"\n")
In [30]:
Out[30]:
The residual should be extremely small (around epsilon machine)
You will perform a benchmark to obtain the assymptotic complexity of your algorithm with respect to the number of unknows in the system.
You will run your algorithm to solve linear systems for matrices of different sizes; you will time the execution time for each of those solves, and you will plot the runtime versus the size of the matrices in a log-log scale. From the plot you will claim the assymptotic complexity of your algorithm with respect to the number of unknowns on the linear
Q2.a You will run the following script, to bechmark the assymptotic complexity of the Julia built-in linear system solver (which is an interface to LAPACK, for more information see https://en.wikipedia.org/wiki/LAPACK and http://www.netlib.org/lapack/)
In [14]:
nSamples = 10;
times = zeros(nSamples,1)
sizes = 2*2.^(0:nSamples-1)
for i = 1:nSamples
m = sizes[i]
# creating the Matrix
A = rand(m,m) + m*eye(m)
# creating the rhs
b = rand(size(A)[1],1)
tic();
x =A\b
times[i] = toc();
end
You will plot the timing using the following script (you will use the layer object to plot two different data-set in the same graph)
In [15]:
using Gadfly
plot(
layer(x = sizes, y = times, Geom.point, Geom.line),
layer(x = sizes, y = 0.00000001*sizes.^3, Geom.point, Geom.line, Theme(default_color=color("red"))),
Scale.y_log10, Scale.x_log10,
Guide.ylabel("Runtime [s]"), # label for y-axis
Guide.xlabel("Size of the Matrix"), # label for x-axis
)
In red you have the cubic scaling and in light blue the numerical runtimes.
What should you expect?
What are you seeing instead?
How can you explain this?
What would happen if you increse the size of the matrices?
Answer: We should expect the two lines to be parallel; however, the zig-zag pattern look like a quadratic complexity. This is a clear example of a prea-assymptotic regime, given that the problem is not big enough we are not able to "see" the leading exponent. If the size of the problems to solve were bigger, we would be able to see the assymptotic scaling.
Q2.b You will modify the script in the question above in order to time the algorithm you wrote.
In [16]:
nSamples = 10;
times2 = zeros(nSamples,1)
sizes2 = 2*2.^(0:nSamples-1)
for i = 1:nSamples
m = sizes2[i]
# creating the Matrix
A = rand(m,m) + m*eye(m)
# creating the rhs
b = rand(size(A)[1],1)
tic();
x = solveGauss(A,b)
times2[i] = toc();
end
In [17]:
using Gadfly
plot(
layer(x = sizes2, y = times2, Geom.point, Geom.line),
layer(x = sizes2, y = 0.00000001*sizes2.^3, Geom.point, Geom.line, Theme(default_color=color("red"))),
Scale.y_log10, Scale.x_log10,
Guide.ylabel("Runtime [s]"), # label for y-axis
Guide.xlabel("Size of the Matrix"), # label for x-axis
)
Based on the runtime scaling you obtained, what is the assymptotic complexity of your function solveGauss?
We can see the two lines are parallel (at least for the last 3 or 4 points) then we can claim that the complexity of the algorithm is cubic as expected.
In [ ]: